C.................... MONPRN.FOR........................................
C....... This program reads data created by the FLIP model calculates
C........ densities production rates etc. and prints them.
      SUBROUTINE MONPRN
      use FIELD_LINE_APEX_PARAMS  !.. altmaX,glatX,glonX,altX,hrX,SIGNLAT

      IMPLICIT DOUBLE PRECISION(A-H,N,O-Z)
      INTEGER ION,NNI,GPRINT,VD,AUR_FILE
      INTEGER NUMUNIT,UNITNUMS(22)     !.. Unit numbers for writing rate warning
      PARAMETER (VD=181)    !.. VD = dimension of vertical grid
      REAL SEC,DEC,ETRAN,F107,F107A,AP,BLON,D,T,CHI,SATN,SATF,SAT,EUV
     > ,UVFAC,GLON,GLONN,GLONF,CHIN,CHIF,GLATN,GLATF,GLATJ,GRD,ZWR
     > ,SIG(2,20),PROB(3,6,37),ZLAM,ZFLUX,SIGABS,SIGION,TPOT,RATFAC
     > ,SW(25),SATX,RE,PLAT,PLON,PI,ANGVEL,DLAT,DLON
     > ,TIMEON,TIMEOF,AURALT,CHIM,SATM,GLONM,GLATM,SIGOXS(22),SIGEXT
     > ,EE,SIGEX(3),SIGI(3),SIGO1D,UTHRS,FLIP_ALT,AE_OP,AE_OX,AE_N2
     > ,AE_O2,AE_NO,AE_TE,AE_TI,GLATIT,RZWR,TSTART,TSTOP,TSAVE,TNEXT
      REAL PEBSC(3),SIGEL(3),HE10830XS,UHED(2)
      REAL EFAC,EFACN,EFACS,ECDT,ECHROM,ECORON   !.. Solar eclipse factors and time step

      DIMENSION D(19),T(2),FD(9),COSDIP(401),HE(401),NT(2),IPRIN(55),
     >          VCON(VD),UE(2,401), Q(401),COOL(22)
      COMMON/RFAC/RATFAC(99)
      COMMON/RK/VC(401)
      COMMON/SIGS/ZFLUX(37),SIGABS(3,37),ZLAM(37),SIGION(3,37),
     > TPOT(3,10),NNI(3),LAMAX
      !-- common to hold the total ion production rates
      REAL SUMION,SUMEXC
      COMMON/BPRTOT/SUMION(3,12,401),SUMEXC(3,12,401)
      COMMON/FJS/N(4,401),TI(3,401),F(20)
      COMMON/MSIS/AP(7),DEC,ETRAN,BLON,F107,F107A,IDAY,SEC,GL(401)
      COMMON/VN/U(2,401),BG(401),BM(401),GR(2,401),R(2,401),SL(401)
      COMMON/ND/ON(401),HN(401),N2N(401),O2N(401),PHION(401),TN(401)
      COMMON/ALT/Z(401),DT,DH,THF,EPS,TF,ITF,JMAX,JMAX1,ITER,ION
      COMMON/SAV/NSAVE(2,401),TISAV(3,401),FY(2,401),UN(401),EHT(3,401)
      COMMON/LPS/SZA(401),EPSN,DC,ISPC,I1,I2,IPV,JTIMAX,IPP,ISKP
      COMMON/ND1/NR(2,VD),ZR(VD),TZ(VD),GZ(VD),BZ(VD),VDT,TE(VD),DZ,JTER
      COMMON/HEPRM/HFY(2,401),HVEL(2,401),HDEN(2,401)
      COMMON/STAW/STR(12,401),RCON(401)
      COMMON/SOL/UVFAC(59),EUV
      COMMON/GPAR/RE,PLAT,PLON,PI,ANGVEL
      COMMON/MIB/OP2D(401),OP2P(401),N2D(401),N4S(401),NNO(401)
      COMMON/VODDN/VN2D(401),VN4S(401),VNNO(401),VO1D(401),IVERT
      COMMON/MINOR/XIONN(9,401),XIONV(9,401),XMASS(9)
      COMMON/VAGRID/JBV,JVB,JNB
      COMMON/CONFAC/IWIND,IVIBN2,IODDN,INPLUS,IHEP,ITEMP,IHPOP
      !-- common used in CMSIS to set longitude for neuts
      COMMON/MJLON/JVERT,JVERC

      DATA TINC , ITSAV,  IFAST , JPRIN, JJJ, ZTRAN, OPLS, DTIMS
     >    /   900,-10000 ,   1   ,  0   ,  0 , 0.0  , 0.0 , 0.0  /
       DATA JWR, JB , ZLO , ZHI ,SETCHI, I556
     >   /  0  ,  1 , 200 , 250 ,   1.7,    0/
      DATA IPRIN/3*1,52*0/,VCON/VD*1/, SW/25*1.0/
      DATA D/19*0.0/, FEXISTS/-9/
      DATA EFAC,EFACN,EFACS,ECDT/4*1.0/   !.. Eclipse factor

      !.. Unit numbers for writing reaction rate warning
      DATA NUMUNIT/7/,UNITNUMS/3,4,6,8,9,10,11,15*0/ 

      READ(5,*) IFAST
      READ(5,*) TSTART
      READ(5,*) TSTOP
      READ(5,*) TINC
      IF(TINC.LT.0.1) TINC=0.1
      TINC=TINC*60
      READ(5,*) ZLO
        IF(ZLO.LT.80) ZLO=80
      READ(5,*) ZHI
        IF(ZHI.LT.ZLO) ZHI=ZLO
      READ(5,*) ZWR
      READ(5,*) IPP
      READ(5,*) IEUVXS
      READ(5,*) ALTP
      DO 171 K=1,49
 171  READ(5,*) IPRIN(K)

      !.. need to print all heating and cooling rates for heating 
      !.. efficiency calculation
      IF(IPRIN(35).NE.0) THEN
        DO I=29,34
          IPRIN(I)=1
        ENDDO
      ENDIF

      DO 213 I=1,55
 213  JPRIN=JPRIN+IPRIN(I)

      CALL RATCHK(NUMUNIT,UNITNUMS) !... check reaction rates standard

      !.. IF TSTART=TSTOP read thru file to find nearest time
      IF(IFAST.NE.1) THEN
         WRITE(6,*) ' @.. Searching for START time on data file'
         CALL FNDTIM(TSTART,TSTOP,TSAVE,IDAY,PCO,BLON)
      ENDIF

      TNEXT=TSTART  !... next time to print

       !.. Use PCO as L-value to determine magnetic latitude and 
      !.. set up Apex Coordinates
      CALL SETUP_APEX(IDAY)

      CALL TFILE(14,AUR_FILE)  !... see if auroral file exists

      JTER=1

      !.. solar declination and ephemeris transit
      CALL SOLDEC(IDAY,SEC/3600.,DEC,ETRAN)

      !..WRITE(8,56)
      WRITE(6,56)
 56   FORMAT(//2X,' ***** CAUTION : check validity of local equilibrium'
     > ,' densities.',/' most are doubtful above 400-500km, NO+ and O2+'
     > ,' are doubtful below 200km at night **'/)

      CALL READ_Q(SCALK, Q)  !.. get coordinate system

      IFLAG=0
C;;;;;;;;;;;;;;;;;;;;;;;; main time step begins ;;;;;;;;;;;;;;;;;;;;;;;;
      DO 221 JTI=1,9999999
      SECSAV=SEC

      !.. Calculate solar eclipse factor in North and South
      CALL ECFUN(SEC,IDAY,45.0,ECDT,ECHROM,ECORON)
      EFACN=ECHROM  !.. 2017-09-08: Use chromospheric emission for photoionization
      IF(EFACN.LT.1) WRITE(6,'(7X,A,F6.3)') 
     >  'Eclipse in Northern Hemisphere. EUV attenuation factor =',EFACN
      CALL ECFUN(SEC,IDAY,-45.0,ECDT,ECHROM,ECORON)
      EFACS=ECHROM  !.. 2017-09-08: Use chromospheric emission for photoionization
      IF(EFACS.LT.1) WRITE(6,'(7X,A,F6.3)') 
     >  'Eclipse in Southern Hemisphere. EUV attenuation factor =',EFACS

      !.. input initial densities and temperatures .......
      CALL RF8DAT(IVLPS,JQ,IWR,PCO,Z0,SCAL,GRD,GLONN,GLONF,GLATN,GLATF,
     >  ZLBDY,FPAS,UV,HFAC,JMAX,IDAY,I1,I2,F107,F107A,AP,BLON,EUV,Z,
     >  UVFAC,RATFAC,TIMEON,TIMEOF,AURALT,PLAT,PLON)
      JMAX1=JMAX-1

      !.. Get EUVAC or SEE EUV fluxes. If UVFAC(58) >=0 use values  
      !.. read from input file. These are user supplied EUV factors
      IF(NINT(UVFAC(58)).EQ.-1.OR.NINT(UVFAC(58)).EQ.-3)
     >   CALL FACEUV(F107,F107A,UVFAC,ZFLUX)
      IF(NINT(UVFAC(58)).EQ.-2.OR.NINT(UVFAC(58)).EQ.-4)
     >   CALL FACSEE(F107,F107A,UVFAC,ZFLUX)

      !.. set up and output the magnetic field line spatial grid .....
       GPRINT=7
       IF(ABS(SEC/3600.0-TSTART).LE.0.01) GPRINT=7
       CALL RFIELD(GPRINT,PCO,RE,Z0,SCAL,SCALK,ZLBDY,Q,VPERP,VTOT
     >           ,COSDIP,JB,JR,JNQ,JTI,1.5D3,ZLO,JVERT,JVERC)

      !-- Find altitude points for setting lat, long for vertical grid
       JALN=1
       JLOB=1
       JHIB=JQ
       JP=JQ
       DO 245 IL=2,JMAX/2
         IF(Z(IL).LT.250) JALN=IL
         IF(ZLO.LT.(Z(IL+1)+Z(IL))/2.AND.ZLO.GE.(Z(IL)+Z(IL-1))/2)
     >      JLOB=IL
         IF(ZHI.LT.(Z(IL+1)+Z(IL))/2.AND.ZHI.GE.(Z(IL)+Z(IL-1))/2)
     >      JHIB=IL
         IF(ALTP.LT.(Z(IL+1)+Z(IL))/2.AND.ALTP.GE.(Z(IL)+Z(IL-1))/2)
     >      JP=IL
 245   CONTINUE

      JALC=JMAX+1-JALN
      !.. see if ALTP closer to next grid point and check limits
      IF(DABS(ALTP-Z(JP+1)).LT.DABS(ALTP-Z(JP))) JP=JP+1
      IF(ALTP.GE.Z(JQ)) JP=JQ
      IF(ALTP.LE.Z(1)) JP=2

      JS=JMAX+1-JP      !.. conjugate point to JP
      !-- get lat and long at 250 km altitude
      CALL AMBS(JALN,Z(JALN),GL(JALN),D,T,CHIN,SATN,GLONN,GLATN)
      CALL AMBS(JALC,Z(JALC),GL(JALC),D,T,CHIF,SATF,GLONF,GLATF)

      !.......... print ............................
      IF(JTI .EQ. 1) THEN
        DO I=3,11
          IF(I.NE.5.AND.I.NE.7) 
     >      CALL PHEADR(I,PCO,GLONN,GLONF,GLATN,GLATF,Z(JALN),Z0
     >        ,JMAX,SCAL,IVIBN2+IODDN)
        ENDDO

        !-- Print the run control parameters 
        CALL RUNPRN(6,ISKP,Z0,HPEQ,FPAS,HFAC,UVFAC,AUR_FILE,PCO)
      ENDIF

       CALL VERGRD(IVLPS,GRD)      !.. set up the vertical grid

      CALL RF8D2(JTX,JK,JC,NMF2N,NMF2F,NOPEQ,NHPEQ,TIEQ,TEEQ,UNN,UNC,
     >  IVLPS,IFLAG,JQ,SEC,N,JMAX,TI,U,STR,RCON,AP,NHEPEQ,UE,F107
     >  ,NNPEQ)

      !... Make sure vertical grid upper boundary is below field line apex
      JBP1=JBV
      DO J=1,JBV
        IF(ZR(J).LT.Z(JQ)) JBP1=J
      ENDDO

      !.. If 3 hour Ap then tell MSIS
      IF(SW(9).GE.0.0.AND.AP(7).GT.0) THEN
        SW(9)=-1.0
        CALL TSELEC(SW)
      ENDIF

      IF(IFLAG.EQ.1) GO TO 233

      DT=SEC-SECSAV
      UTHRS=SEC/3600.0
      IF(JTI.EQ.1) DT=0.0001

      !.. Call AMBS to get local times SATN, SATF and CHI
      CALL AMBS(JALN,Z(JALN),GL(JALN),D,T,CHIN,SATN,GLONN,GLATN)
      CALL AMBS(JALC,Z(JALC),GL(JALC),D,T,CHIF,SATF,GLONF,GLATF)

      IF(I556.NE.-9) WRITE(8,556)
      IF(I556.NE.-9) WRITE(9,556)
      I556=-9
 556  FORMAT(/1X,'** This file gives the quantities along the magnetic'
     > ,1X,'field line grid.' 
     > ,/1X,'O+, H+, He+, N+ are from full calculations. IF [He+] and'
     > ,1X,'[N+] = 0 they were not calculated'
     > ,/1X,'If odd N diffusion equations are solved they are'
     > ,1X,'interpolated onto the field line grid.'
     > ,/1X,'Otherwise N2D and N4S are from local equilibrium. N(4S) is'
     > ,1X,'from MSIS'
     > //'**** NO OUTPUT? try the files for vertical altitude grid ***')

      KKK=KKK+1
      IF(MOD(KKK,11).EQ.0.OR.KKK.EQ.1) WRITE(6,317)
      IF(IFAST.NE.0.AND.KKK.EQ.1) WRITE(3,317)

 317     FORMAT(/10X,'[------Northern hemisphere-------]',5X,
     >   '[------Southern hemisphere------]',/3X,'  UT    LT  SZA'
     >   ,3X,'NmF2   HmF2  N2* Lat Lon     LT  SZA   NmF2   HmF2'
     >   ,2X,'N2* Lat Lon Kp  F107  JTI DelT L-Sh YYYYddd')
 117     FORMAT(F8.2,F6.2,I4,1P,E9.2,I5,0P,F5.2,2I4,2X,F6.2,I4,1P
     >   ,E9.2,I5,0P,F5.2,2I4,I3,I5,I6,I5,F6.3,I8)
      RVLOC=RCON(JK)
      IF(RVLOC.LT.1) RVLOC=(RCON(JK-1)+RCON(JK+1))*0.5
      IF(RVLOC.LT.1) RVLOC=1.0
      RVCON=RCON(JC)
      IF(RVCON.LT.1) RVCON=(RCON(JC-1)+RCON(JC+1))*0.5
      IF(RVCON.LT.1) RVCON=1.0
      KPIND=0
      IF(AP(2).GT.2) KPIND=INT(ALOG(AP(2))*1.7-0.84)
      !.. get F2 peak electron density but only if a true O+ peak. NMF2N and
      !.. NMF2F are already read in but not the hmF2s. 1st Northern hemisphere
      HMF2N=Z(JK)
      NMF2N=N(1,JK)+N(3,JK)
      IF((N(1,JK-1)+N(3,JK-1) .LT. N(1,JK)+N(3,JK)).AND.
     >     (N(1,JK+1)+N(3,JK+1) .LT. N(1,JK)+N(3,JK))) 
     > CALL F2PEAK(N(1,JK-1)+N(3,JK-1),N(1,JK)+N(3,JK)
     >      ,N(1,JK+1)+N(3,JK+1),Z(JK-1),Z(JK),Z(JK+1),NMF2N,HMF2N)
      !.. Southern hemisphere
      HMF2F=Z(JC)
      NMF2F=N(1,JC)+N(3,JC)
      IF((N(1,JC-1)+N(3,JC-1) .LT. N(1,JC)+N(3,JC)).AND.
     >     (N(1,JC+1)+N(3,JC+1) .LT. N(1,JC)+N(3,JC))) 
     > CALL F2PEAK(N(1,JC-1)+N(3,JC-1),N(1,JC)+N(3,JC)
     >      ,N(1,JC+1)+N(3,JC+1),Z(JC-1),Z(JC),Z(JC+1),NMF2F,HMF2F)

       CALL ACTUAL_DAY(IDAY,SEC,JDAY,SX)    !.. gets actual day of year

       WRITE(6,117) UTHRS,SATN,NINT(57.3*CHIN),NMF2N
     >   ,NINT(HMF2N),RVLOC,NINT(GLATN),NINT(GLONN),SATF
     >   ,NINT(57.3*CHIF),NMF2F,NINT(HMF2F),RVCON,NINT(GLATF)
     >   ,NINT(GLONF),KPIND,NINT(F107),JTI,NINT(DT),PCO,JDAY

       IF(IFAST.NE.0) WRITE(3,117) UTHRS,SATN,NINT(57.3*CHIN),NMF2N
     >   ,NINT(HMF2N),RVLOC,NINT(GLATN),NINT(GLONN),SATF
     >   ,NINT(57.3*CHIF),NMF2F,NINT(HMF2F),RVCON,NINT(GLATF)
     >   ,NINT(GLONF),KPIND,NINT(F107),JTI,NINT(DT),PCO,JDAY

      !... Check to see whether to skip printing
      IF(NMF2N.LE.0.OR.NMF2F.LE.0.OR.TIEQ.LE.0.OR.TEEQ.LE.0) GO TO 221
      IF(IFAST.NE.0) GO TO 221
      IF(UTHRS.LT.TSTART-0.006) GO TO 221

      !.. check to see if we need to stop
      IF(UTHRS.GT.TSTOP+.0049) GO TO 233
      IF(ABS(TSTOP-TSTART).LT.0.01.AND.UTHRS.GT.TSTOP) GO TO 233

      !.. test to see if enough time has passed for the next print
      IF(UTHRS.LT.TNEXT-TINC/360000.) GO TO 221
      TNEXT=TNEXT+TINC/3600.0
      JJJ=JJJ+1

      N(1,JQ)=NOPEQ
      N(2,JQ)=NHPEQ
      N(4,JQ)=NHEPEQ
      XIONN(1,JQ)=NNPEQ
      TI(2,JQ)=TIEQ
      TI(3,JQ)=TEEQ

      DO 5 I=1,2
 5    NT(I)=0.0

      DO 734 J=1,JMAX
      VC(J)=RCON(J)
      TI(1,J)=TI(2,J)
      !..... total tube content ........
      IF(J.GT.1) THEN
         NT(1)=NT(1)+BM(JB)*N(1,J)*(SL(J)-SL(J-1))/BM(J)
         NT(2)=NT(2)+BM(JB)*N(2,J)*(SL(J)-SL(J-1))/BM(J)
      ENDIF
 734   CONTINUE


      JFQ=JMAX-JNQ+1
      FYNOP=N(1,JNQ)*U(1,JNQ)*BM(JB)/BM(JNQ)
      FYNHP=N(2,JNQ)*U(2,JNQ)*BM(JB)/BM(JNQ)
      FYFOP=N(1,JFQ)*U(1,JFQ)*BM(JB)/BM(JFQ)
      FYFHP=N(2,JFQ)*U(2,JFQ)*BM(JB)/BM(JFQ)

      !.. prints headers in vertical files
      CALL PRINT_VER_HEAD(IVLPS,JJJ,JPRIN,Z(JNQ),IPRIN(1))

      CALL AMBS(JK,Z(JK),GL(JK),D,T,CHIM,SATM,GLONM,GLATM) !.. get MSIS quantities

      WRITE(3,356) UTHRS,SATM,NINT(57.3*CHIM),NINT(GLATM),NINT(GLONM)
     > ,NINT(HMF2N),NINT(UNN/100.),NMF2N,FYNHP,FYNOP,NINT(TI(2,JK))
     > ,NINT(TI(3,JK)),NINT(T(2)),D(2),D(4),D(3)

      !... call AMBS for neutral densities
      CALL AMBS(JC,Z(JC),GL(JC),D,T,CHIM,SATM,GLONM,GLATM) !.. get MSIS quantities

      WRITE(4,356) UTHRS,SATM,NINT(57.3*CHIM),NINT(GLATM),NINT(GLONM)
     > ,NINT(HMF2F),(NINT(-UNC/100.)),NMF2F,-FYFHP,-FYFOP,NINT(TI(2,JC))
     > ,NINT(TI(3,JC)),NINT(T(2)),D(2),D(4),D(3)
 356  FORMAT(F8.3,F6.2,5I5,1P,E9.2,2E10.2,3I6,1P,9E9.2)

      !.. prints headers in field line files
      CALL PRINT_FL_HEAD(JJJ,JPRIN,Z(JP),Z(JS))

      !-- This section for writing field line densities and production rates.
      !.. northern hemisphere ....
      FYOPJP=N(1,JP)*U(1,JP)*BM(1)/BM(JP)
      FYHPJP=N(2,JP)*U(2,JP)*BM(1)/BM(JP)
      HFLUX1=7.7E+5*TI(3,JP)**2.5*(TI(3,JP+1)-TI(3,JP-1))/
     > (SL(JP+1)-SL(JP-1))
      HIFLX1=1.15E+4*TI(2,JP)**2.5*(TI(2,JP+1)-TI(2,JP-1))/
     > (SL(JP+1)-SL(JP-1))

      CALL AMBS(JP,Z(JP),GL(JP),D,T,CHI,SATX,GLON,GLATJ) !.. get MSIS quantities
      !.. added 2023-01-09 to output Hedin winds UHED(1),UHED(2)
      CALL hwm14(JDAY,SX,REAL(Z(JP)),GLATJ,GLON,SATX,F107A,F107,AP,UHED)

      WRITE(9,361) UTHRS,SATX,NINT(57.3*CHI),N(1,JP),N(2,JP)
     > ,NINT(TI(2,JP)),NINT(TI(3,JP)),FYOPJP,FYHPJP,D(2),D(4),D(3)
     > ,D(7),D(1),NINT(T(2)),HFLUX1,N(4,JP)
     > ,JDAY,KPIND,NINT(F107),NINT(F107A),UNN/100.0,UHED(1),UHED(2)
       
      !.. print North major and minor neutral densities at an altitude as a function of time
      CALL TVPRINT(1,JMAX/2,JBP1,Z(JP),GRD,GL(JVERT),ISOL,MAXIT,IVLPS)

      !.. southern hemisphere ....
      FYOPJS=-N(1,JS)*U(1,JS)*BM(1)/BM(JS)
      FYHPJS=-N(2,JS)*U(2,JS)*BM(1)/BM(JS)
      HFLUX2=-7.7E+5*TI(3,JS)**2.5*(TI(3,JS+1)-TI(3,JS-1))/
     > (SL(JS+1)-SL(JS-1))
      HIFLX2=-1.15E+4*TI(2,JS)**2.5*(TI(2,JS+1)-TI(2,JS-1))/
     > (SL(JS+1)-SL(JS-1))

      CALL AMBS(JS,Z(JS),GL(JS),D,T,CHI,SATX,GLON,GLATJ) !.. get MSIS quantities
      !.. added 2023-01-09 to output Hedin winds UHED(1),UHED(2)
      CALL hwm14(JDAY,SX,REAL(Z(JS)),GLATJ,GLON,SATX,F107A,F107,AP,UHED)

      WRITE(8,361) UTHRS,SATX,NINT(57.3*CHI),N(1,JS),N(2,JS)
     > ,NINT(TI(2,JS)),NINT(TI(3,JS)),FYOPJS,FYHPJS,D(2),D(4),D(3)
     > ,D(7),D(1),NINT(T(2)),HFLUX2,N(4,JS)
     > ,JDAY,KPIND,NINT(F107),NINT(F107A),UNC/100.0,UHED(1),UHED(2)

      !.. print South major and minor neutral densities at an altitude as a function of time
      CALL TVPRINT((JMAX+2)/2,JMAX,JBP1,Z(JP),GRD,GL(JVERC),ISOL,MAXIT,
     > IVLPS)

      IF(IPRIN(23).NE.0) GO TO 177
      IF(JPRIN-IPRIN(1)-IPRIN(37).LE.0.AND.ZWR.EQ.0) GO TO 371

 177  CONTINUE
 
      JPEAK=(JMAX+1)/2
      DO 17 J=1,JMAX
      !..... obtain neutral densities -- ambs calls msis ......
      CALL AMBS(J,Z(J),GL(J),D,T,CHI,SAT,GLON,GLATJ)
      HE(J)=D(1)
      ON(J)=D(2)
      HN(J)=D(7)
      N2N(J)=D(3)
      O2N(J)=D(4)
      TN(J)=T(2)
      SZA(J)=CHI
      N4S(J)=D(8)
 17   CONTINUE

      !------------------- end neutral atmosphere loop -------------------

      !-- transfer vertical grid odd N to field line grid
      IF(JPRIN.GT.0.OR.IPRIN(37).GT.0) CALL TRSTR(IVLPS)

      !.. This section allows FLIP to substitute measured values
      !.. from a file for FLIP O+, O, and N2. The file data are
      !.. interpolated onto the FLIP grid. If the data are bad,
      !.. zero values are returned and FLIP values are retained.
      !.. O, O2, and N2 are replaced in AMBS
      DO J=1,JMAX
         FLIP_ALT=Z(J)
         GLATIT=GLATN
         IF(J.GT.JMAX/2) GLATIT=GLATF
         CALL GET_AE_DENS(FLIP_ALT,UTHRS,GLATIT
     >      ,AE_OP,AE_OX,AE_N2,AE_O2,AE_NO,AE_TE,AE_TI)
         IF(AE_OP.GT.100) N(1,J)=AE_OP
         IF(AE_NO.GT.10) VNNO(J)=AE_NO
         IF(AE_TI.GT.300) TI(1,J)=AE_TI
         IF(AE_TI.GT.300) TI(2,J)=AE_TI
         IF(AE_TE.GT.300) TI(3,J)=AE_TE
      ENDDO

      !.. determine h+ flux by simple method
      IF(IPRIN(55).NE.0.AND.PCO.GT.1.2) THEN
         IF(ZTRAN.LE.0.0) THEN
            CALL PHEADR(111,PCO,GLONN,GLONF,GLATN,GLATF,Z(JALN),Z0
     >           ,JMAX,SCAL,IVIBN2+IODDN)
           !.. WRITE(111,251)
         ENDIF
      CALL FLXLIMIT(JR,JQ,JK,Z,N,TI,BM,GR,ON,HN,TN,SL,JZ,FYLIM,HPLOSS,
     >   SH_OPLUS,ZTRAN)
         WRITE(111,256) UTHRS,SATN,NINT(57.3*CHIN),FYNHP,FYLIM,HPLOSS
     >   ,NINT(ZTRAN)
      ENDIF
 251  FORMAT(/4X,'UT',4X,'LT',2X,'Chi',3X,'PhiFLIP',3X,'PhiLIM',3X,
     >   'ILOSS     Crit')
 256  FORMAT(F7.2,F6.2,I4,1P,3E10.2,I7)
      !.. Determine HmF2 by diffusion=loss
      IHMF2=0
      IF(IHMF2.NE.0) THEN
         DIP=ACOS(COSDIP(JK))
         OXJ=1.4E-7*F107/70.
         CHID=CHIN
         UM=UNN/100.
         IF(OPLS.LE.0) OPLS=N(1,JK)
         OSAV=OPLS
         DTIM=SEC-DTIMS
         DTIMS=SEC
         CALL HMF2(I,ZS,Z(JK),ON(JK),N2N(JK),O2N(JK),TN(JK),TI(3,JK)
     >      ,UM,DIP,CHID,OXJ,OPLS,OSAV,DTIM)
         WRITE(11,356) UTHRS,SATN,NINT(57.3*CHIN),NINT(Z(JK)),NINT(ZS)
     >     ,NINT(UNN/100.),NMF2N,OPLS,FYNHP,FYLIM,HPLOSS,DTIM,OSAV,DIP
         WRITE(11,356) UTHRS,SATN,NINT(57.3*CHIN),NINT(Z(JK)),NINT(ZS)
     >      ,NINT(UNN/100.),NMF2N,ON(JK),N2N(JK),O2N(JK),TN(JK)
     >      ,TI(3,JK),UM,DIP,CHID
      ENDIF

      IF(JPRIN-IPRIN(1)-IPRIN(37).LE.0.AND.ZWR.EQ.0) GO TO 485

      !... Obtain symmetric ambient conditions
      IF(IVLPS.EQ.-1.OR.IVLPS.EQ.-2) CALL SYM(JMAX,HE,EHT,UN,Z,SZA)

      !..... evaluate primary production rates from primpr
      DO J=1,JMAX
        IF(J.LE.JMAX/2) EFAC=EFACN
        IF(J.GT.JMAX/2) EFAC=EFACS
        CALL PRIMPR(J,Z,ON,N2N,O2N,HE,SZA,TN,JMAX,SETCHI,IVLPS,EFAC)
      ENDDO

      !..  pe2s=2-stream prog to get electron heating rate;
       IF(CHIN.LE.3.0.OR.CHIF.LE.3.0) CALL PE2S
     >   (Z,ON,N2N,O2N,HE,BM,BG,JMAX,DH,N,TI,EHT,TN,SZA,FPAS,2,ZWR,1,SL)

      !..... If there is an EUV flux file need to reevaluate primary rates
      IF(FEXISTS.EQ.-9) CALL TFILE(42,FEXISTS)
      IF(FEXISTS.NE.0) THEN
        DO J=1,JMAX
          IF(J.LE.JMAX/2) EFAC=EFACN
          IF(J.GT.JMAX/2) EFAC=EFACS
          CALL PRIMPR(J,Z,ON,N2N,O2N,HE,SZA,TN,JMAX,SETCHI,IVLPS,EFAC)
        ENDDO
      ENDIF

      !.. AURSUB calculates auroral electron production rates
      RZWR=ZWR
      CALL AURSUB(1,JMAX/2,UTHRS,Z,ON,O2N,N2N,N,TI,EHT,DTMAX,RZWR)

      !-- Sum the EUV, photoelectron, and auroral production rates
      CALL SUMPRD(JMAX)

      !.. call routine to print heating factors ....
      IF(ZWR.LT.0) CALL HTPRIN(JJJ,JPRIN,9,UTHRS,SATN,CHI,HFLUX1,HFLUX2
     > ,HIFLX1,HIFLX2)

      DO 27 J=1,JMAX
      !.. 1st term due to unknown heat source. HFAC=1 -> HF=1E9 eVcm-3s-1
      IF(Z(J).GT.3000) EHT(3,J)=3.0*HFAC/PCO**4 + EHT(3,J)
      !N(3,J)=0.0
      IF(Z(J).LT.Z0.OR.Z(J).GT.900) GO TO 27
      !.. call cminor to get no+, o+(2d), n2+ & o+(2p) densities,
      CALL CMINOR(0,J,0,IPR,FD,9,HE(J))
      !-- PHION=total O+(4S) prod, including EUV, e*, dissoc of O2 and minor ions
      PHION(J)=SUMION(1,7,J)+SUMION(2,4,J)+SUMION(2,5,J)+FD(9)
      N(3,J)=XIONN(1,J)+FD(1)+FD(2)+FD(3)+FD(4)+FD(5)
 27   CONTINUE

      !... Obtain symmetric ambient conditions
      IF(IVLPS.EQ.-1.OR.IVLPS.EQ.-2) CALL SYM(JMAX,HE,EHT,UN,Z,SZA)

      IF(IVIBN2+IODDN.GT.0) THEN
        !.. Printing odd N, Minor ions etc.-- Northern hemisphere
        IF (JPRIN.GT.0) CALL VPRINT(JALN,JPRIN,IPRIN,IVLPS,1,JMAX/2
     >   ,GL(JVERT),ISOL,MAXIT,JBP1,ZLO,ZHI,IPP,Z(JP),GRD)

      !-- Southern hemisphere
        IF (JPRIN.GT.0) CALL VPRINT(JALC,JPRIN,IPRIN,IVLPS,(JMAX+2)/2,
     >   JMAX,GL(JVERC),ISOL,MAXIT,JBP1,ZLO,ZHI,IPP,Z(JP),GRD)
      ENDIF
      !.............................................................

 485  CONTINUE
      DO 385 IS=2,36
        JP=1
      DO 385 J=1,JMAX/2+1
        IF(J.LT.JLOB.OR.J.GT.JHIB) GO TO 385
        IF(MOD(J+IPP-1,IPP).NE.0.AND.J.NE.1) GO TO 385
        !..IF(Z(J).LT.100) GO TO 385
        IF(IPRIN(IS).EQ.1) CALL CMINOR(IS,J,JP,0,FD,9,HE(J))
        JU=JMAX+1-J
        IF(IPRIN(IS).EQ.1) CALL CMINOR(IS,JU,JP,0,FD,8,HE(JU))
        JP=0
 385  CONTINUE
 371  CONTINUE

      !.. print densities and temperatures from subroutine loops....
      IF(IPRIN(1).EQ.0) GO TO 877
      IF(IPRIN(1).NE.0) WRITE(9,859)
      IF(IPRIN(1).NE.0) WRITE(8,859)
 859  FORMAT(/'  J    ALT   TE   TI    [O+]     [H+]      VO+'
     >  ,8X,'VH+     [e]     SZA   [He+]      VHe+      N+'
     >  ,7X,'VN+     O+ +N2v     BMag        SL'
     >  ,8X,'W_Bmer    HWMmer  HWMzon')
      DO 854 J=1,JQ
      JU=JMAX+1-J
      IF(J.LT.JLOB.OR.J.GT.JHIB) GO TO 854
      IF(J.EQ.1.OR.J.EQ.JQ) GO TO 852
      IF(MOD(J+IPP-1,IPP).EQ.0) GO TO 852
      GO TO 854
 852  CONTINUE

      NEDEN=N(1,J)+N(2,J)+N(3,J)+N(4,J)+XIONN(1,J)
      CALL AMBS(J,Z(J),GL(J),D,T,CHI,SAT,GLON,GLATJ)
      !.. added 2023-01-09 to output Hedin winds UHED(1),UHED(2)
      CALL hwm14(JDAY,SX,REAL(Z(J)),GLATJ,GLON,SAT,F107A,F107,AP,UHED)
      WRITE(9,888) J,NINT(Z(J)),NINT(TI(3,J)),NINT(TI(1,J)),N(1,J),
     > N(2,J),U(1,J),U(2,J),NEDEN,CHI*57.3,N(4,J),UE(1,J)
     > ,XIONN(1,J),XIONV(1,J),VC(J),BM(J),SL(J)
     > ,-UN(J)/COSDIP(J)/100.0,UHED(1),UHED(2)   !.. added 2023-01-09
      SEDEN=N(1,JU)+N(2,JU)+N(3,JU)+N(4,JU)+XIONN(1,JU)
      CALL AMBS(JU,Z(JU),GL(JU),D,T,CHI,SAT,GLON,GLATJ)
      !.. added 2023-01-09 to output Hedin winds UHED(1),UHED(2)
      CALL hwm14(JDAY,SX,REAL(Z(JU)),GLATJ,GLON,SAT,F107A,F107,AP,UHED)
      WRITE(8,888) JU,NINT(Z(JU)),NINT(TI(3,JU)),NINT(TI(1,JU)),N(1,JU)
     > ,N(2,JU),U(1,JU),U(2,JU),SEDEN,CHI*57.3,N(4,JU)
     > ,UE(1,JU),XIONN(1,JU),XIONV(1,JU),VC(JU),BM(JU),SL(JU)
     > ,-UN(JU)/COSDIP(JU)/100.0,UHED(1),UHED(2)  !.. added 2023-01-09 
 854  CONTINUE
 888  FORMAT(I4,I6,I6,I5,1P,2E9.2,1P,2E10.2,1P,E9.2,0P,F6.1,1P,E9.2
     > ,1P,E10.2,1P,E9.2,1P,E10.2,1P,2E10.2,1P,E13.4,0P,3F9.2
     > ,1P,22E10.2)

 877  IF(IPRIN(29).NE.1) GO TO 221
      DO 458 I=1,2
      DO 458 J=1,JQ
      JU=J
      IF(I.EQ.2) JU=JMAX+1-J
      IF(J.LT.JLOB.OR.J.GT.JHIB) GO TO 458
      IF(J.EQ.1.OR.J.EQ.JQ) GO TO 258
      IF(MOD(J+IPP-1,IPP).EQ.0) GO TO 258
      GO TO 458
 258  CONTINUE
      CALL ECRATS(6+I,JU,Z(JU),TI(3,JU),TI(2,JU),TN(JU),N(1,JU),N(2,JU)
     > ,N(3,JU),ON(JU),N2N(JU),O2N(JU),HN(JU),HE(JU),EHT(3,JU),TLSS
     > ,FOL,COOL)
 458  CONTINUE

 221  CONTINUE
 233  CONTINUE

      !.. print EUV fluxes and cross sections
      IF(IEUVXS.NE.1) GO TO 819
         CALL PHEADR(17,PCO,GLONN,GLONF,GLATN,GLATF,Z(JALN),Z0
     >       ,JMAX,SCAL,IVIBN2+IODDN)
         CALL PARAMS(1,LMAX)
         CALL PROBS(1,PROB,ZLAM,LMAX,NNI)
      WRITE(17,790)
 790  FORMAT(/22X,'electron impact cross sections for N2'
     > ,/3X,'E     A3      B3pi      C3      a1pi    a1sig ',
     > '    b1pi    b1sig      W3     B3sig  a''1sigu  w1del')
      EE=0.5
      DELTE=1
      DO JE=1,1000
        EE=EE+DELTE
        IF(EE.GT.99) DELTE=10
        IF(EE.GT.999) GO TO 801
        CALL N2FOB(EE,SIG,IKS)
        WRITE(17,101) EE,(SIG(1,IK),IK=3,5),(SIG(1,IK),IK=7,IKS)
      ENDDO
 801  CONTINUE

      WRITE(17,792)
 792  FORMAT(/9X,'Total electron impact cross sections (excitation,'
     > ,1X,'ionization, elastic) AND atomic oxygen partial cross'
     > ,1X,'sections'
     > ,/3X,'E     O*       O2*      N2*      O+      O2+      N2+'
     > ,5X,'O_el      O2_el   N2_el'
     > '     O(1D)      O(1S)     1304      989      1027    1356',
     > '    8446    7774   He10830')
      EE=0.5
      DELTE=1
      DO JE=1,1000
        EE=EE+DELTE
        IF(EE.GT.99) DELTE=10
        IF(EE.GT.999) GO TO 819
        CALL OXSIGS(EE,SIGOXS,SIGEXT)
        CALL SIGEXS(EE,SIGEX,SIGI,SIGO1D)
        CALL ELASTC(EE,PEBSC,SIGEL)
        WRITE(17,101) EE,(SIGEX(IK),IK=1,3),(SIGI(IK),IK=1,3)
     >  ,(SIGEL(IK),IK=1,3)
     >  ,(SIGOXS(IK),IK=1,6),SIGOXS(3)*0.6,SIGOXS(6)*0.5,HE10830XS(EE)
      ENDDO
 819  CONTINUE

      CALL FLIP_MODS   !.. Print FLIP updates

      CLOSE (UNIT=3)
      STOP

 101  FORMAT(F5.1,1P,22E9.2)
 361  FORMAT(F8.3,F6.2,I4,1P,2E9.2,I5,I6,1P,2E10.2,1P,5E9.2,I5
     > ,1P,E9.1,1P,E8.1,I8,3I5,0P,3F9.2)
      END
C:::::::::::::::::::::::: PRINT_FL_HEAD :::::::::::::::::::::
C..... Prints hemisphere specific headers in Field Line grid files
C..... Modified 2023-01-09 for Hedin winds UHED(1),UHED(2) output
      SUBROUTINE PRINT_FL_HEAD(JJJ,JPRIN,ALTN,ALTS)
      DOUBLE PRECISION ALTN,ALTS
      IF(JJJ.EQ.1) WRITE(9,363) ALTN
 363  FORMAT(/4X,'NORTH ION DENSITIES, TEMPS, HWM14 Winds AT ALT='
     >   ,F9.0,' KM, W_Bmer = Wind @ hmF2 along magnetic meridian')
      IF(JJJ.EQ.1) WRITE(8,364) ALTS
 364  FORMAT(/4X,'NORTH ION DENSITIES, TEMPS, HWM14 Winds AT ALT='
     >   ,F9.0,' KM, W_Bmer = Wind @ hmF2 along magnetic meridian')
      IF(JJJ.EQ.1.OR.JPRIN.GT.0) WRITE(9,365)
      IF(JJJ.EQ.1.OR.JPRIN.GT.0) WRITE(8,365)
 365  FORMAT(/4X,'UT',5X,'LT',2X,'CHI',3X,'[O+]',5X,'[H+]',5X,'Ti   Te '
     > ,4X,'FO+',7X,'FH+',6X,'[O]',6X,'[O2]',5X,'[N2]',5X,'[H]',6X
     > ,'[He]',5X,'Tn',3X,'HFLX   [He+]    IDAY     KP F107 F107A'
     > ,2X,'W_Bmer  HWMmer   HWMzon')
      RETURN
      END
C:::::::::::::::::::::::: PRINT_VER_HEAD :::::::::::::::::::::
C... Prints hemisphere specific headers in vertical grid files
      SUBROUTINE PRINT_VER_HEAD(IVLPS,JJJ,JPRIN,ALT,IPRIN1)
      DOUBLE PRECISION ALT

      IF(JJJ.NE.1.AND.IVLPS.LE.0) RETURN

      IF(JJJ.EQ.1.AND.IVLPS.LE.0) WRITE(3,470)
      IF(JJJ.EQ.1.AND.IVLPS.LE.0) WRITE(4,470)
 470    FORMAT(//2X,' *?*?* There are no vertical profiles in this file'
     >  ,1X,'because there are no odd N densities on the FLIP run data'
     >  ,1X,'file'/)

      IF(JJJ.EQ.1) WRITE(3,351) NINT(ALT)
      IF(JJJ.EQ.1) WRITE(4,352) NINT(ALT)

 351  FORMAT(/6X,'NORTHern hemisphere parameters at hmF2'
     > ,1X,'(WIND + northward),'
     > ,1X,'and plasmaspheric refilling fluxes (Phi** @',I5,' km)')
 352  FORMAT(/6X,'SOUTHern hemisphere parameters at hmF2'
     > ,1X,'(WIND + northward),'
     > ,1X,'and plasmaspheric refilling fluxes  (Phi** @',I5,' km)')

      IF(JJJ.EQ.1.OR.JPRIN-IPRIN1.GT.0) WRITE(3,451)
      IF(JJJ.EQ.1.OR.JPRIN-IPRIN1.GT.0) WRITE(4,451)
 451  FORMAT(/4X,'UT    LT  Chi  GLAT GLON HmF2 Wind   NmF2     PhiH+'
     > ,5X,'PhiO+     Ti    Te    Tn     [O]     [O2]     [N2]')
      RETURN
      END
